/** Author: David Powell
This file produces Figure 4 + Figure 6B + Figure A7 + Figure A8 + Figure A10 + Figure A13b + Figure A17 + Figure A18 + Figure A19 + Figure A20 + Figure C1 + Figure C2 + Figure C3
It also includes the results for Table A2
**/


clear all
set more off
set mat 800
set seed 8721

global dir "/jules/b/dpowell"
global DATA "${dir}/purdue/replication/DATA"
global OUTPUT "${dir}/purdue/replication/output"


use ${DATA}/finaldata

gen overdose_rate=100000*(overdose/totpop)
gen opioid_rate=100000*opioids/totpop

gen overdose_no_cocaine_rate=100000*(overdoses_no_cocaine/totpop)
gen opioid_no_cocaine_rate=100000*opioids_no_cocaine/totpop
gen cocaine_no_opioids_rate=100000*cocaine_no_opioids/totpop

gen opioid_not406_rate=opioid_rate
replace opioid_not406_rate=100000*opioids_not406/totpop if year>1998



gen heroin_rate=100000*(heroin/totpop)
gen t402_rate=100000*(t402/totpop)
gen t404_rate=100000*(t404/totpop)

gen cld_rate=100000*cld/totpop
gen suicide_rate=100000*(suicide-suicide_overdose)/totpop

gen post=(year>1995&year<2001)+2*(year>2000&year<2011)+3*(year>2010)
tab stfips, gen(ss) nof
tab year, gen(ttt) nof
drop ttt13
foreach nnn of numlist 1983/1994 1996/2017 {
	gen tt_`nnn'=nontripl*(year==`nnn')
}



***TIME SERIES GRAPHS***
preserve
gcollapse *_rate [aw=totpop], by(nontripl year) fast
# delimit ;


twoway (connected overdose_rate year if nontripl==1, color(red) msize(vtiny) lwidth(thick) sort yaxis(1))   (connected overdose_rate year if nontripl==0 , color(blue) lpattern(dash) msize(vtiny) lwidth(thick)  sort yaxis(1)), graphregion(color(white)) xline(1996) xlabel(1983(2)2017, angle(45))   legend(label(1 "Non-Triplicate States") label(2 "Triplicate States")  rows(1) order(1 2)) 
scheme(s2mono)  ylabel(,nogrid)
xtitle("") ytitle("Deaths per 100,000" ,  axis(1));
gr export ${OUTPUT}/fig4a.eps, replace;

twoway (connected opioid_rate year if nontripl==1, color(red) msize(vtiny) lwidth(thick)  sort yaxis(1))   (connected opioid_rate year if nontripl==0, color(blue) lwidth(thick)  lpattern(dash) msize(vtiny) sort yaxis(1)), graphregion(color(white)) xline(1996) xlabel(1983(2)2017, angle(45))   legend(label(1 "Non-Triplicate States") label(2 "Triplicate States")  rows(1) order(1 2)) 
scheme(s2mono)  ylabel(,nogrid)
xtitle("") ytitle("Deaths per 100,000" ,  axis(1));
gr export ${OUTPUT}/fig4c.eps, replace;


# delimit cr
restore

local aaa=1/7
local nnn=1
***FIGURE 4 + Figure A17 (A and B) + Figure A18
foreach var of varlist overdose_rate opioid_rate suicide_rate cld_rate overdose_no_cocaine_rate opioid_no_cocaine_rate cocaine_no_opioids_rate opioid_not406_rate {
	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0

	reg `var' ss* tt_* ttt*  [aw=totpop], cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}
	/* These variables save the Table A2 results */
	gen dbeta`nnn'=0
	gen dlow`nnn'=0
	gen dhigh`nnn'=0
	gen dpv`nnn'=0


	****TABLE A2****
	qui reg `var' ss* tt_1996-tt_2017 ttt*  [aw=totpop] if year>1990, cluster(stf) 
	replace dbeta`nnn'=(_b[tt_1996]+_b[tt_1997]+_b[tt_1998]+_b[tt_1999]+_b[tt_2000])/5 if post==1
	replace dbeta`nnn'=(_b[tt_2001]+_b[tt_2002]+_b[tt_2003]+_b[tt_2004]+_b[tt_2005]+_b[tt_2006]+_b[tt_2007]+_b[tt_2008]+_b[tt_2009]+_b[tt_2010])/10 if post==2
	replace dbeta`nnn'=(_b[tt_2011]+_b[tt_2012]+_b[tt_2013]+_b[tt_2014]+_b[tt_2015]+_b[tt_2016]+_b[tt_2017])/7 if post==3
	boottest .2*tt_1996+.2*tt_1997+.2*tt_1998+.2*tt_1999+.2*tt_2000=0, boottype(wild) weighttype(webb) reps(9999) seed(23857389)
	matrix A=r(CI)
	replace dlow`nnn'=A[1,1] if post==1
	replace dhigh`nnn'=A[1,2] if post==1
	replace dpv`nnn'=r(p) if post==1
	boottest .1*tt_2001+.1*tt_2002+.1*tt_2003+.1*tt_2004+.1*tt_2005+.1*tt_2006+.1*tt_2007+.1*tt_2008+.1*tt_2009+.1*tt_2010=0, boottype(wild) weighttype(webb) reps(9999) seed(23857389)
	matrix A=r(CI)
	replace dlow`nnn'=A[1,1] if post==2
	replace dhigh`nnn'=A[1,2] if post==2
	replace dpv`nnn'=r(p) if post==2
	boottest `aaa'*tt_2011+`aaa'*tt_2012+`aaa'*tt_2013+`aaa'*tt_2014+`aaa'*tt_2015+`aaa'*tt_2016+`aaa'*tt_2017=0, boottype(wild) weighttype(webb) reps(9999) seed(23857389)
	matrix A=r(CI)
	replace dlow`nnn'=A[1,1] if post==3
	replace dhigh`nnn'=A[1,2] if post==3
	replace dpv`nnn'=r(p) if post==3

	boottest (tt_1996+tt_1997+tt_1998+tt_1999+tt_2000=0) (tt_2001+tt_2002+tt_2003+tt_2004+tt_2005+tt_2006+tt_2007+tt_2008+tt_2009+tt_2010=0) (tt_2011+tt_2012+tt_2013+tt_2014+tt_2015+tt_2016+tt_2017=0), boottype(wild) weighttype(webb) reps(9999) seed(23857389)
	local nnn=`nnn'+1
}

preserve
bys year: keep if _n==1
# delimit ;
twoway(rarea low1 high1 year ,  sort color(gs14))  (connected beta1 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1));
gr export ${OUTPUT}/fig4b.eps, replace;


twoway(rarea low2 high2 year,  sort  color(gs14))  (connected beta2 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1));
gr export ${OUTPUT}/fig4d.eps, replace;

twoway(rarea low3 high3 year,  sort  color(gs14))  (connected beta3 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA17a.eps, replace;


twoway(rarea low4 high4 year,  sort  color(gs14))  (connected beta4 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA17b.eps, replace;

twoway(rarea low5 high5 year,  sort  color(gs14))  (connected beta5 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA18a.eps, replace;

twoway(rarea low6 high6 year,  sort  color(gs14))  (connected beta6 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA18b.eps, replace;

twoway(rarea low7 high7 year,  sort  color(gs14))  (connected beta7 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA19.eps, replace;

twoway(rarea low8 high8 year,  sort  color(gs14))  (connected beta8 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figC2.eps, replace;

# delimit cr
restore


****FIGURE A20*******
local nnn=12
*create state-specific trends
foreach var of varlist ss* { 
	gen trend`var'=(year-1982)*`var'
}
foreach var of varlist overdose_rate opioid_rate {
	reg `var' ss* ttt* trend* [aw=totpop] if year<1996
	predict res_`var', res	

	gen beta`nnn'=0
	gen se`nnn'=.	
	gen low`nnn'=0
	gen high`nnn'=0
	qui reg res_`var' ss* tt_* ttt*  [aw=totpop], cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		replace se`nnn'=_se[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}
	local nnn=`nnn'+1
}
drop trend*

preserve
bys year: keep if _n==1
# delimit ;
twoway (connected beta1 year , color(black) msize(vtiny) sort yaxis(1))  (connected beta12 year , color(green) lpattern(dash) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Overdoses, State Trends") label(1 "Overdoses")  rows(1) order(1 2))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA20a.gph, replace);
gr export ${OUTPUT}/figA20a.eps, replace;
twoway (connected beta2 year , color(black) msize(vtiny) sort yaxis(1))  (connected beta13 year , color(green) lpattern(dash) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Opioid Overdoses, State Trends") label(1 "Opioid Overdoses")  rows(1) order(1 2))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA20b.gph, replace);
gr export ${OUTPUT}/figA20b.eps, replace;
# delimit cr
restore

drop beta* high* low*

****FIGURE 6b*****
preserve
replace nontripl=0 if triplold==1
drop tt_*
foreach nnn of numlist 1983/1994 1996/2017 {
	gen tt2_`nnn'=triplold*(year==`nnn')
	gen tt_`nnn'=nontripl*(year==`nnn')
}
local nnn=4
gen beta`nnn'=0
gen b2eta`nnn'=0

	reg overdose_rate tt_* tt2_* ttt* ss*  [aw=totpop], cluster(stf) 
	foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		replace b2eta`nnn'=_b[tt2_`nnn2'] if year==`nnn2'
	}

# delimit ;

twoway (connected beta4 year if tripl==1 & year>1982 & year<2018, color(red)  lwidth(thick) msize(vtiny) sort yaxis(1)) 
 (connected b2eta4 year if tripl==1 & year>1982 & year<2018, color(green)  lwidth(thick) lpattern(shortdash_dot) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) yline(0, lpattern(dash) lcolor(gs6)) xline(1996) xlabel(1983(2)2017, angle(45))  legend(label(2 "Former Triplicates") label(1 "Never Triplicates")  rows(1) order(2 1))  ylabel(,nogrid) xtitle("Year") ytitle("Estimated Coefficient"  ,  axis(1)) ;
gr export ${OUTPUT}/fig6b.eps, replace;



# delimit cr
restore


local nnn=3
****FIGURE A7*******
foreach var of varlist overdose_rate opioid_rate {
	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0
	qui reg `var' ss* tt_* ttt* xx* [aw=totpop], cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}
	local nnn=`nnn'+1
}



foreach var of varlist overdose_rate opioid_rate {
	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0
	qui reg `var' ss* tt_* ttt*, cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}
	local nnn=`nnn'+1
}



foreach var of varlist overdose_rate opioid_rate {
	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0
	qui reg `var' ss* tt_* ttt* xx*, cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}
	local nnn=`nnn'+1
}

drop ttt*
egen regionxtime=group(region year)
tab regionxtime, gen(ttt) nof
foreach var of varlist overdose_rate opioid_rate {
	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0
	qui reg `var' ss* tt_* ttt* xx* [aw=totpop], cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}
	local nnn=`nnn'+1
}
drop ttt*
tab year, gen(ttt) nof

preserve
bys year: keep if _n==1
# delimit ;

twoway(rarea low3 high3 year,   sort color(gs14))  (connected beta3 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA7e.gph, replace);
gr export ${OUTPUT}/figA7e.eps, replace;

twoway(rarea low4 high4 year,   sort color(gs14))  (connected beta4 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA7f.gph, replace);
gr export ${OUTPUT}/figA7f.eps, replace;

twoway(rarea low5 high5 year,   sort color(gs14))  (connected beta5 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA7a.gph, replace);
gr export ${OUTPUT}/figA7a.eps, replace;

twoway(rarea low6 high6 year,  sort  color(gs14))  (connected beta6 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA7b.gph, replace);
gr export ${OUTPUT}/figA7b.eps, replace;

twoway(rarea low7 high7 year,   sort color(gs14))  (connected beta7 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA7c.gph, replace);
gr export ${OUTPUT}/figA7c.eps, replace;

twoway(rarea low8 high8 year,  sort  color(gs14))  (connected beta8 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA7d.gph, replace);
gr export ${OUTPUT}/figA7d.eps, replace;

twoway(rarea low9 high9 year,   sort color(gs14))  (connected beta9 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA8a.gph, replace);
gr export ${OUTPUT}/figA8a.eps, replace;

twoway(rarea low10 high10 year,  sort  color(gs14))  (connected beta10 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA8b.gph, replace);
gr export ${OUTPUT}/figA8b.eps, replace;

# delimit cr
restore


***FIGURE A13b
foreach var of varlist overdose_rate  {
	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0
	qui reg `var' ss* tt_* ttt* [aw=totpop] if nontripl==0 | stf==21 | stf==26 | stf==18 | stf==53 | stf==46, cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}
	local nnn=`nnn'+1
}

preserve
bys year: keep if _n==1
# delimit ;
twoway(rarea low11 high11 year,  sort  color(gs14))  (connected beta11 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) saving(figA13b.gph, replace);
gr export ${OUTPUT}/figA13b.eps, replace;
# delimit cr
restore


****FIGURE A10***
drop beta* low* high*
local nnn=1
foreach var of varlist heroin_rate t402_rate t404_rate {
	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0
	reg `var'  tt_* ttt* if year>1998 [aw=totpop] , cluster(stf) 
	foreach nnn2 of numlist 1999/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}
	local nnn=`nnn'+1
}
preserve
bys year: keep if _n==1
# delimit ;
twoway(rarea low1 high1 year if  year>1998, sort  color(gs14))  (connected beta1 year if  year>1998, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1999(1)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA10b.eps, replace;

twoway(rarea low2 high2 year if year>1998, sort color(gs14))  (connected beta2 year if year>1998, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1999(1)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA10a.eps, replace;

twoway(rarea low3 high3 year if year>1998, sort  color(gs14))  (connected beta3 year if year>1998, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1999(1)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA10c.eps, replace;

# delimit cr
restore


***FIGURE A17 (detrended)
local nnn=1
gen trend=tripl*(year-1983)
drop beta* low* high*
foreach var of varlist suicide_rate cld_rate {

	reg `var' ttt* trend ss*  [aw=totpop] if year<1996
	gen a_`var'=`var'-trend*_b[trend]

	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0

	reg a_`var' ss* tt_* ttt*  [aw=totpop], cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	
	}

	local nnn=`nnn'+1
}
drop trend
preserve
bys year: keep if _n==1
# delimit ;
twoway(rarea low1 high1 year ,  sort color(gs14))  (connected beta1 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))  ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1));
gr export ${OUTPUT}/figA17c.eps, replace;
gr export figA17c.pdf, replace;

twoway(rarea low2 high2 year,  sort  color(gs14))  (connected beta2 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1));
gr export ${OUTPUT}/figA17d.eps, replace;
gr export figA17d.pdf, replace;




# delimit cr
restore





****FIGURE C1****

foreach yyy of numlist 1983/1994 1996/2017 {
	gen zz_`yyy'=NTR*(year==`yyy')
}
drop beta* low* high*
local nnn=1
foreach var of varlist overdose_rate opioid_rate {
	gen beta`nnn'=0
	gen low`nnn'=0
	gen high`nnn'=0
	gen b2eta`nnn'=0
	gen l2ow`nnn'=0
	gen h2igh`nnn'=0
	qui reg `var' ss* tt_* zz* ttt* [aw=totpop], cluster(stf) 
	qui foreach nnn2 of numlist 1983/1994 1996/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'	

		replace b2eta`nnn'=_b[zz_`nnn2'] if year==`nnn2'
		boottest zz_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace l2ow`nnn'=A[1,1] if year==`nnn2'
		replace h2igh`nnn'=A[1,2] if year==`nnn2'
	}
	local nnn=`nnn'+1
}

preserve
bys year: keep if _n==1
# delimit ;
twoway(rarea low1 high1 year,  sort  color(gs14))  (connected beta1 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figC1a.eps, replace;

twoway(rarea l2ow1 h2igh1 year,  sort  color(gs14))  (connected b2eta1 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figC1b.eps, replace;

twoway(rarea low2 high2 year,  sort  color(gs14))  (connected beta2 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figC1c.eps, replace;

twoway(rarea l2ow2 h2igh2 year,  sort  color(gs14))  (connected b2eta2 year , color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1))   ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figC1d.eps, replace;
# delimit cr
restore




***FIGURE C3***
drop post
gen post=(year>1995)*nontripl
levelsof stf, local(sss)
gen beta=.
gen high=.
gen low=.
foreach num of numlist `sss' {
	qui reg overdose_rate ss* post ttt*  [aw=totpop] if stf!=`num', cluster(stf) 
	replace beta=_b[post] if stf==`num'
	boottest post, boottype(wild) weighttype(webb) reps(9999) seed(23857389)
	matrix A=r(CI)
	replace low=A[1,1] if stf==`num'
	replace high=A[1,2] if stf==`num'
}
egen statenum=group(stf)

label def stf 6 California 4 Arizona 41 Oregon 32 Nevada 1 Alabama 5 Arkansas 8 Colorado 2 Alaska 9 Connecticut 10 Delaware 11 D.C. 12 Florida 13 Georgia 15 Hawaii 16 Idaho 17 Illinois 18 Indiana 19 Iowa 20 Kansas 21 Kentucky 22 Louisiana 23 Maine 24 Maryland 25 Massachusetts 26 Michigan 27 Minnesota 28 Mississippi 29 Missouri 30 Montana 31 Nebraska 33 "New Hampshire" 34 "New Jersey" 35 "New Mexico" 36 "New York" 37 "North Carolina" 38 "North Dakota" 39 Ohio 40 Oklahoma 41 Oregon 42 Pennsylvania 44 "Rhode Island" 45 "South Carolina" 46 "South Dakota" 47 Tennessee 48 Texas 49 Utah 50 Vermont 51 Virginia 53 Washington 54 "West Virginia" 55 Wisconsin 56 Wyoming
label val stf stf
# delimit ;
twoway (scatter beta statenum) (rcap high low statenum), saving(leaveoneout.gph, replace)  ylabel(,nogrid) graphregion(color(white))  ytitle("Estimate") xtitle("Excluded State") legend(label(1 "Estimated Coefficients") label(2 "95% Confidence Interval")  rows(1) order(1 2))
xlabel(1 "Alabama" 2 "Alaska" 3 "Arizona" 4 "Arkansas" 5 "California" 6 "Colorado"
7 "Connecticut"
8 "Delaware"
9 "District of Columbia"
10 "Florida"
11 "Georgia"
12 "Hawaii"
13 "Idaho"
14 "Illinois"
15 "Indiana"
16 "Iowa"
17 "Kansas"
18 "Kentucky"
19 "Louisiana"
20 "Maine"
21 "Maryland"
22 "Massachusetts"
23 "Michigan"
24 "Minnesota"
25 "Mississippi"
26 "Missouri"
27 "Montana"
28 "Nebraska"
29 "Nevada"
30 "New Hampshire"
31 "New Jersey"
32 "New Mexico"
33 "New York"
34 "North Carolina"
35 "North Dakota"
36 "Ohio"
37 "Oklahoma"
38 "Oregon"
39 "Pennsylvania"
40 "Rhode Island"
41 "South Carolina"
42 "South Dakota"
43 "Tennessee"
44 "Texas"
45 "Utah"
46 "Vermont"
47 "Virginia"
48 "Washington"
49 "West Virginia"
50 "Wisconsin"
51 "Wyoming", angle(90) labsize(vsmall));
graph export ${OUTPUT}/figC3.eps, replace;
# delimit cr












